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Abstract 

Non-invasive marks including pigmentation patterns, acquired scars, and 
genetic markers, are often used to identify individuals in mark-recapture exper- 
iments. Such data present new statistical challenges including that the marks 
may change over time, that the marks may be misidentified at non-negligible 
rates, and that only part of the population may possess a specific mark. Moti- 
vated by the analysis of data from the ECOCEAN online whale shark catalog, 
we consider issues that arise when individuals may be identified from multi- 
ple, non-invasive marks. Past approaches to this problem include modeling 
data from only one mark and combining inferences obtained from each mark 
separately under the assumption that marks are independent. We describe a 
Bayesian method that makes use of the data from multiple marks and properly 
accounts for the dependence between marks from the same individual. Results 
from the analysis of the ECOCEAN data and simulation studies are included 
and show that the new method provides more precise inference than analyzing 
the data from one mark while maintaining nominal coverage rates for credible 
intervals. 

Keywords: Latent multinomial model; Mark-recapture; Multiple marks; Non-invasive 
marks; Photo-identification; Whale sharks 



20 1 Introduction 



21 Non-invasive marks (also called natural marks) include patterns in pigmentation, ge- 

22 netic markers, acquired scars, or other natural characertistics that allow researchers 

23 to identify individuals in a population without physical capture. Such marks have 

24 been used for a long time in mark-recapture studies of some classes of animals, par- 

25 ticularly marine mammals, and are now being used more widely. The advantages of 

26 non-invasive marks over man-made marks are first that non-invasive marks can be ob- 

27 served passively from a distance or through the collection of secondary material (hair 

28 samples or scat) so that individuals do not have to be handled and second, depending 

29 on the type of mark, that every individual is marked from birth. However, non- 
30 invasive marks also present several modeling challenges including that marks may be 

31 misidentified at non-negligible rates, that individuals' marks may change over time, 

32 and that certain types of marks may be restricted to a subset of the population. We 

33 consider the problem of modeling about the demographics of a population based on 

34 mark-recapture data when individuals may be identified from multiple non-invasive 

35 marks. 

36 The specific application we consider concerns the aggregation of whale sharks in 

37 Ningaloo Marine Park (NMP), off the west coast of Australia, each year between April 

38 and July. During this time, individual whale sharks are located by tour companies and 

39 photographs are taken by tourists and tour operators who upload their images to the 

40 online ECOCEAN whale shark library. Whale sharks can be identified by the unique 

41 pattern of spots on their flanks, and computer assisted methods are used to match 

42 photographs in the library and generate capture histories which provide information 

43 about the timing of the sharks' arrival and departure from NMP and their survival 



44 across years (Holmberg et al. 2009). Yoshizaki et al. (2009) and Yoshizaki et al. 



45 (2011) list many more applications involving the use of non-invasive marks with wide 



range of species including photographs of large cats (cheetahs, snow leopards, and 
tigers), scar patterns on marine mammals (manatees and whales), skin patterns of 
reptiles and amphibians (snakes, crocodiles, and salamanders), and genetic marks 
(bears, wombats, and whales). 

Though non-invasive marks provide advantages they also present modeling chal- 



51 lenges that do not arise with man-made tags. Da-Silva et al. (2003) and Da-Silva 



52 (2006) consider difficulties that occur if only a portion of the population possesses 



the non-invasive mark; for example, if individuals are identified from body scars that 
they acquire over time. Some individuals are effectively uncatchable because they 
cannot be identified, and abundance estimates must be inflated to account for this 
portion of the population. A related challenge is that non-invasive marks may change 



jpver time, as individuals acquire new scars or their skin patterns change. |YoshizakT 



et al. (2011 ) calls these as "evolving marks" and develops a least squares based method 



that allows for estimation of population demographics under the assumption that it is 
impossible to match the mark before and after the change creating a "ghost capture 



6i history" ( Yoshizaki et al. 2011 pg. 29). 



Work has also been done on the problem of misidentification of non-invasive marks, 



63 particularly for genetic marks. Lukacs and Burnham (2005) develops an extension of 



the Jolly-Seber model to account for misidentification of genetic markers under the 
assumption that the same error cannot occur twice and that errors do not produces 
spurious matches with other individuals. This produces an apparent excess of indi- 



67 _viduals captured once from which the rate of misidentification is estimated. Yoshizaki 



et al. (2011) shows that the estimators of Lukacs and Burnham (2005) are biased 



69 even in large samples, and presents an alternative least squares approach based on 



to the same assumptions. Further, Wright et al. (2009) presents a method that relaxes 



71 the assumptions that errors in identification when the probability that markers are 

72 misread can be estimated from multiple tests of each genetic sample. 

73 The problem we address arises when a single individual can be identified from 

74 more than one non-invasive mark, but the marks are not always linked. This occurs 

75 in our example because a single whale shark may photographed from either the right 

76 or the left side, but the spot patterns on each side are distinct and cannot be linked 

77 without further information. This means that it may not be possible to connect ob- 

78 servations of a shark photographed one from the left and later from the right and to 

79 know that the same individual was if fact encountered twice. Naively modeling the 
so observed encounter histories will artificially inflate the number of individuals observed 
si and create dependence between the encounter histories, violating a key assumption of 
82 most mark-recapture models. One solution is to construct encounter histories from a 



single mark, Holmberg et al. (2009) considers only left-side photographs, but this may 



84 remove a large number of histories from the data. Alternatively, Wilson et al. ( 1999 ) 

85 combines inferences from each mark post hoc by averaging separates estimates of abun- 
se dance and computing standard errors under the assumption that these estimates are 

87 independent. It is simple to show that the final point estimate has similar properties 

88 to the individual estimates - the overall estimate will be unbiased/consistent if each of 

89 the separate estimates is unbiased/consistent - but the assumption of independence 



90 is violated and the resulting standard errors will be biased low. Madon et al. (2011) 

91 describes another method to estimate abundance by adjusting the sufficient statistics 

92 required to compute the Jolly-Seber estimator. While the approach seems plausible, 

93 we have concerns with the implementation; though the observed counts underesti- 



94 mate some of the statistics and overestimate others, Madon et al. (2011) uses a single 



95 adjustment factor and constrains its value to be between and 1. Bonner (2012) 



96 discusses these issues and shows that the correct adjustment factor is not constant 



97 and may be greater than 1 for some statistics. Moreover, the simulations |Madon et al. 



98 (2011) presents indicate a clear problem in that the coverage of confidence intervals 

99 is much lower than their nominal value, even when the population is large and the 

100 capture probability is close to 1. 

101 Our objectives are to construct a model that accounts for the dependence between 

102 multiple non-invasive marks and provides valid inference about the population's de- 

103 mographics while making use of all available data. We do so by constructing an 

104 explicit model of the observation process that allows for multiple marks to be ob- 

105 served on different occasions and then applying Bayesian methods of inference via 
loe Markov chain Monte Carlo (MCM) sampling. Our model can be cast as a special 



107 case of the latent multinomial model (LMM) of Link et al. (2010). Section [2] describes 

108 the whale shark data and Section [3] gives details on the structure of our model and 

109 our methods of inference. We provide results from a simulation study and the appli- 
no cation of our method to the whale shark data in Sections [4] and |5j and conclude with 
in a discussion of the method and further extensions in Section [6] 



n2 2 Data 

in Data for our analysis was obtained from the ECOCEAN on-line whale shark library 

114 (available at |www. w hale shark. org[ ). This library contains photographs of whale 

us sharks taken by recreational divers and tour operators worldwide and submitted on- 

n6 line. Computer assisted algorithms are then used to identify the pattern of spots 

in on the shark's flank and scan the database for matches with previous photographs 

us which are confirmed visually. Quality control studies indicated that the probability 



n9 of a mismatch is approximately .02 (Holmberg et al. , 2009). The library has been 



120 operational since 2003, and over 41,000 photographs have been submitted by more 
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121 than 3300 contributors. Further details on the study site, the observation protocols, 



122 and the algorithms for matching photographs are provided in Holmberg et al. (2009). 

123 We model only the data collected from the northern ecotourism zone of NMP 

124 during the 16 week period between April 1 and July 31, 2008, which we have divided 

125 into 8 capture periods of 2 weeks each. Sharks may be encountered multiple times 

126 during a single capture period and so five possible events may occur. In each 2 week 

127 period, a shark may: 

128 1) not be encountered at all (event 0) 

129 2) be photographed from the left only (event L), 

130 3) be photographed from the right only (event R), 

131 4) be photographed from both sides simultaneously on at least on encounter (event 

132 S), or 

133 5) be photographed from both sides but never simultaneously (event B). 

134 Problems with identification arise because it is only possible to match the skin pat- 

135 terns from the two sides of a shark if photographs of both sides are taken simultane- 

136 ously during at least one capture period - i.e., there is at least one S in its encounter 

137 history. Otherwise, an individual that was photographed from both sides will con- 

138 tribute two encounter histories to the data set - one containing the observations of 

139 its right side and the other containing the observations of its left side. 
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3 Methods 



141 3.1 Modeling 

142 The challenge in modeling this data is that some encounter histories are unobservable 

143 so that the recorded histories may not reflect the true histories of the encountered 

144 individuals. Suppose that an individual's true encounter history is 00L0B0R0. This 

145 history is not observable because the two sides of the individual were never pho- 

146 tographed simultaneously and the marks on its right and left sides cannot be linked. 
14? Instead, the individual will contribute two observed histories to the data - 00L0L000 

148 and 0000R0R0. Working in the opposite direction, if the observed data includes the 

149 histories 00L0L000 and 0000R0R0 then it is not possible to know if these histories 

150 come from one individual encountered on three occasions or from two individuals 

151 encountered twice each. 

152 With K — 8 capture occasions there are 5 8 — 1 = 390, 624 possible true encounter 

153 histories (we ignore the zero history because our model conditions on an individual 

154 being captured once at least). Of these, 325,558 (83.3%) can occur in the observed 

155 data because they contain photographs from the left-side only, photographs from the 
we right-side only, or at least one simultaneous encounter. The remaining 65,026 (16.7%) 
is? contain photographs taken from both the left and right sides, including via event B, 

158 but never simultaneously and are unobservable. In practice, we need not consider 

159 all 65,026 unobservable histories. Given a set of observed histories, the possible 
wo true histories can be constructed by combining each history containing only left- 
lei side photographs with each history containing only right-side photographs as follows: 

162 when event L appears in the first history and appears in the second enter L in the 

163 combined history, when appears in the first history and R in the second enter R in 

164 the combined history, and when L appears in the first and R in the second then B is 
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165 entered into the combined history. 

lee To account for the uncertainty in the true encounter histories when modeling the 



167 demographics of the population we adapt the LMM of Link et al. (2010). This model 

168 applies when the latent counts of the true histories follow a multinomial distribution 
leg and determine the observed counts through a known, linear relationship. Suppose 
170 that L unique histories were observed, let rij be the number of times that the j th 
m history was observed, and set n = (m, . . . , ni) T . Further, suppose that there are 

172 L* true histories that could have generated the observed histories, let n*j be the 

173 number of individuals with the j th true history, and set n* = (n*, . . . ,n* L *) T . The 

174 LMM supposes first that n = An* for some fixed matrix A and second that n* has 

175 a multinomial distribution whose density, f(n*\n*,0), depends on n*, representing 

176 either the total population size or the number of unique individuals encountered, 

177 and the vector of parameters 0. The likelihood for this model is proportional to the 

178 probability of n given both n* and 6 and is computed by summing f(n*\n*, 6) over 

179 all possible configurations of n* that satisfy equation (??) and produce the correct 

180 value of n* (n* = J2 n *j)- Evaluating this sum directly is generally not possible, and, 



lei instead, Link et al. (2010) applies Bayesian inference treating n* as a random variable 



182 and sampling from the joint posterior distribution 

n(n* ,n*,6\n) oc l(n = An*)f(n*\n*, 6)n(n*, 6) (1) 

183 where 7r(n*, 6) is the chosen prior for n* and 0. 

184 To apply the LMM to our problem, we need to identify the set of possible true 

185 histories that could have generated the observed data and the relationship between 
lee the corresponding counts. Suppose that the observed data contains L\ unique histo- 
187 ries including only left-side photographs, L 2 unique histories including only right-side 
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188 photographs, and L 3 unique histories with at least one encounter of both sides simul- 

189 taneously so that L = Li + L 2 + L 3 . Let Wi be the LixK matrix whose rows contains 

190 those histories including only left-side photographs and rii the Li-vector of observed 

191 counts for each row of W\. Define W 2 , n 2 and W 3 , n 3 similarly for those histories 

192 including only right side photographs and those including at least one simultaneous 

193 encounter so that the complete vector of observed counts is n = (rii T , n 2 T , n 3 T ) T . 

194 The additional entries in the set of possible true histories are then generated by com- 

195 bining each row of W\ with each row of W 2 , as described above. Let W 4 denote 
we the L 4 x K matrix of unobservable histories (L 4 = LiL 2 ) whose Li(i — 1) + j row is 
19? formed by combining row i of W\ with row j of W 2 . Let denote the number of 

198 individuals encountered whose true encounter history given by the Li{i — 1) + j th row 

199 of W 4 and n 4 = (n 411 , . . . ,nl Lll , . . . ,nl 1L2 , . . . ,nl LlL2 ) T so that the complete vector 

200 of true counts is n* = (n^ T , n 2 T , nf- ', n 4 T ) T '. 

201 The set of linear constraints mapping n* to n is then obtained by considering 

202 how each true history contributes to each observed history. Specifically, the number 

203 of times that a history in W\ or W 2 , say w, is observed must equal the sum of the 

204 counts for all of the true histories that could have generated w. These constraints 

205 are determined by the set of L linear equations 




i=i 

Li 




(2) 



n 3k = n* 3kj , 



k = l,...,L 3 



206 one for each observed history. Note that there are L equations and L* unknowns, 

207 implying that there are L* — L — L\ free parameters. It is most convenient to choose 
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208 n\ as the set of free parameters, and the possible values can be bounded by noting 

209 that < n^j < min(nij, n%j). This limits the possible configurations of n\ and makes 



210 it easier to sample from the posterior distribution, as described in Section |3.3 

2n As an example, suppose that the observed data comprises the six histories listed 

212 in the top of Table [I] The possible true histories that generated this data include 

213 the six observed histories plus the four unobservable histories generated by combining 

214 the histories including only of left-side photographs (histories 1&2) with the histo- 

215 ries including only right-side photographs (histories 3&4), as given in the bottom of 

216 the table. The relationship between n and n* is then determined by the six linear 
21? equations 

n\ = n\ + n* 7 + n* 9) n 3 = + n* 7 + n%, n 5 = n\ 
n 2 = n* 2 + n* 8 + n* , n 4 = n\ + n* 9 + n* l0l n 6 = n%. 

218 The free parameters and these values are bounded so that < rij < 

219 min(ni,ns), < nl < min(n2,Ti3), < nl < min(ni,n4), and < n\ < min(n2,n4). 

220 [Table 1 about here.] 

221 3.2 Model of Observed Encounter Histories 

222 This framework for modeling data from multiple non-invasive marks is general and 

223 can incorporate any of the standard mark-recapture models for the true encounter 

224 histories. We are most interested in the timing of the sharks' arrival and departure 

225 from NMP and have developed an extension of the Link-Barker- Jolly-Seber (LBJS) 



226 model (Link and Barker, 2005) which accounts for the presence of multiple marks. 



227 The LBJS model describes the demographics of a population through two sets of 



10 



228 parameters: 1) the recruitment rate, k = 1, . . . , K — 1, defined as the number of 

229 individuals that enter the population between occasions k and k + 1 per individual 

230 present on occasion k and 2) the survival probability, <pk, k — 1, . . . , K — 1, defined as 

231 the probability that an individual present on occasion k is also present on occasion 

232 k + 1. For the whale sharks, (pk corresponds to the probability that an individual 

233 remains at NMP rather than true survival. It is assumed that all emigration is 

234 permanent, that individuals behave independently, and that the apparent survival 

235 probability does not depend on how long an individual has been present (or any 

236 other factors). 

237 The process of observing individuals is then modeled in terms of two further sets 

238 of parameters: 1) the probability that an individual is captured/encountered given 

239 that it is present in the study area and 2) the conditional probabilities of each event. 

240 Specifically, we define pt to be the probability that an individual is photographed 

241 at least one time during occasion k given that it is present, k = 1, . . . , K, and pj 

242 to be the probability that the individual is photographed from the left-side only 

243 (j = 1), from the right-side only (j = 2), from both sides simultaneously at least once 

244 (j = 3), or from both sides but never simultaneously (j = 4) during period k. This 

245 last set of parameters is what distinguishes our model of the true histories from the 



246 original model of Link and Barker (2005). Our model assumes that all individuals 

247 present on an occasion have the same probability of capture, that encounters are 

248 independent between individuals and across the occasions, and that the conditional 

249 event probabilities do not vary between individuals or over time. 
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3.3 Inference 



251 As noted by Link et al. (2010), the LMM lends itself to Bayesian inference via MCMC 

252 sampling with data augmentation, and the Bayesian specification of the model is com- 

253 pleted by assigning prior distributions to all of the parameters. We have selected non- 
254 informative, hierarchical prior distributions where possible. Specifically, the survival 

255 and capture probabilities were assigned prior distributions 

logit(0 fc ) ~jV(^,<rJ), k = l,...,K-l 
logit(p fc ) ~ N(fi p ,al), k = l,...,K 

256 with 

/i^, fM p ~ N(0, 2) and 0^, o p ~ HT(3, .9) 

257 where HT denotes the half t-distribution with 3 degrees of freedom and scale pa- 

258 rameter .9. These hyperparameters were chosen because they provide distributions 

259 that are close to uniform on (0, 1) for both the hierarchical means, \x$ and fi p , and 

260 the base demographic parameters, <pk and p^- The recruitment probabilities were 

261 assigned prior distributions of the form 

log(/ fc )~iV(0/,*/), k = l,...,K-l 

262 with 

H f ~ N(0, .25) and a f ~ HT{3, .9), 

263 and the number of unique individuals encountered was assigned the discrete uni- 

264 form prior on 0,...,M for some M guaranteed not to be smaller than n*. This 

265 is satisfied by the data dependent prior n* ~ U{0, . . . , L x + L 2 + L 3 }, because 
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266 n* = Ylj=i n j — X^'=i n j' but the exact value does not affect our computations. 

267 Finally, the conditional event probabilities were assigned a Dirichlet prior with pa- 

268 rameter a — I4 implying that p±, p 2 , P3, and pi were marginally distributed a priori 

269 with identical beta distributions having mean .25 and probability .9 between .00 and 

270 .86. All prior distributions were assumed independent. 

271 An MCMC algorithm for sampling from the distribution in ([!]) can then be im- 

272 plemented as follows. Letting n*( k \ n*( k \ and 6^ k \ denote the values sampled on the 

273 k th iteration, the values on the next iteration are generated in two steps: 

274 1. Update n*( fc ) and n*^ given 0^ by: 

275 a) setting n*( curr ) = n*^ and 7r*( curr ) = n*^ k \ 

276 b) updating n* A11 by: 

277 i) setting n*' = n* (curr ) 

278 ii) sampling a new value for n\ n from {0, . . . , min(nn, 7221)} according to some 

279 proposal distribution, g(n4 n |n*( curr) , 6^), 

280 iii) ensuring the constraints in equation ^ remain satisfied by setting 

*' \ *' 

Li 

"21 = n 2 - J^nii, and 



8=1 



n In . 



iv) rejecting n*' if n*^ < or n* 2 \ < and otherwise accepting n*' and setting 
n *(curr) _ n *' anc | n *(curr) _ n *' w ^ probability: 



a(n*, n* ) = min < 1, 



f(n*' 


\n*' } 0( k ))q(n*£ mi) \ 


n*',6^) 


j( n *(curr)| 


n *(curr );6> (fe)) g ( n *' ii 


| n *(curr) ? #(*:)) 
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283 c) repeating step b) for each n*^-, i — 1, . . . , Li, j — 1, . . . , L 2 , and finally 



284 d) setting n*( fc+1 ' = n *( curr ) and n*( fc+1 ) = n H curr J. 

285 2. Update 6^ given n*( fc+1 ) and n*( fc+1 ) through a series of Metropolis-Hastings (MH) 

286 steps appropriate for the selected model. 

287 The steps for updating each element of n* in this algorithm have an intuitive inter- 

288 pretation. In step i), a new count is proposed for one of the unobservable histories. 

289 In step ii), the counts of the corresponding observed histories in Wi and W2 are ad- 

290 justed to maintain the correct observed counts. The new proposal is rejected immedi- 

291 ately if this produces negative counts, and otherwise, a MH accept-reject step is con- 

292 ducted. One choice for q{n\^\n* , 6) is the discrete uniform on {0, ... , min(nij, n 2 j)} 

293 in which case the proposal, n*' is accepted whenever the likelihood is increased, 

294 f{n*'\n*',0) > /(n< curr )|n*( curr ),6>). 

295 This MCMC sampling algorithm was implemented in JAGS using the rjags in- 



[ Plummer 


2003, 


2011 


Team 


2012 



297 were run in parallel in order that convergence could be monitored with the Gelman- 



298 Rubin-Brooks (GRB) diagnostic (Brooks and German, 1998) as implemented in the 



299 R package CODA (Plummer et al. , 2006). One challenge we faced was obtaining 



initial values for different sets of parameters were dispersed but consistent with one 
other. To do this, we generated initial values for each chain in three steps. First, we 
simulated random values of <pk and fk, k = 1, . . . , K — 1, under one of three schemes - 
1) logit(0 fe ) ~ iV(logit(.9),.l) and log(/) ~ iV(log(.l), .1), 2) logit(0 fc ) ~ N(, .1) and 
log(/) ~ iV(log(.5),.9), 3) logit(0 fe ) ~ iV(logit(.l),.l) and log(/) ~ iV(log(.9), .1). 
We then obtained values of pk, k = 1,...,K, by fitting the LBJS model to the 
data obtained from the left-side photographs only but treating the sampled values 
of 4>k and fk as fixed. Finally, we obtained values of n* by fitting the our model 
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308 to the full data set but treating all of and pk as fixed. An R package which 

309 implements these methods and provides functions to format the data and to call 

310 JAGS through the rjags interface is available from the website of the first author at 

311 |www. simon. bonners . ca/MultiMark, 



312 4 Simulation Study 

313 To assess the performance of the model presented in the previous section we conducted 

314 simulation studies under a variety of scenarios. Here we present the results from three 

315 simulation scenarios which illustrate our main results. 

316 In the first scenario, we compared the performance of the new model (henceforth 

317 the two-sided model) with two alternative models. First, we considered models using 

318 reduced data sets containing only data from the left or right side photographs (the 

319 one-sided models). Capture histories were constructed by combining all events that 

320 include a left-side photograph, namely L, S, and B, and all right-side photographs 

321 were ignored, or vice versa. The models we fit to this data were equivalent to the 



322 LBJS model with prior distributions as given in Section 3.3 Second, we developed a 



323 Bayesian equivalent of the method of combining inferences from the two sides under 



324 the assumption of independence as in (Wilson et al. 1999) (combined inference). To 



325 do this, we fit separate models to the data from the left and right side photographs and 

326 then averaged the values sampled on each iteration of the separate MCMC samplers 

327 prior to computing summary statistics. For example, let (f)^'^ and represent 

328 the values of 4>t drawn on the k th iterations of the MCMC samplers run separately for 

329 models of the the left and right-side data and Var^(^) and Var^ (0 t ) be the separate 

330 posterior variances. Combined inference for (f> t was obtained by first computing the 
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331 inverse variance weighted average of the samples drawn on the k th iteration 

.(*) = (Var^(0 t )0f' L) +VarW(0 f )#' fi) ) 
Var( L )(0 t )+Var( R )(0 t ) 

332 and then computing summary statistics from the new chain 0^,0^, . . .. The ad- 

333 vantage of this method is that posterior standard deviations and credible intervals 

334 can be computed directly from the new chain without distributional assumptions. 

335 Note that the mean of the values in the new chain is exactly equal to the inverse 

336 variance weighted average of means from the separate chains. We expected that, in 

337 general, posterior standard deviations and credible intervals from the one-sided mod- 

338 els would be larger/wider than the corresponding standard deviations and intervals 

339 from the two-sided model, and that posterior standard deviations and credible inter- 

340 vals produced by combined inference would be smaller /narrower than those from the 

341 two-sided model but would not achieve the nominal coverage probability. 

342 One hundred data sets for this scenario were generated under the assumption that 

343 each event was equally likely conditional on capture (pi = P2 = P3 = P4 = -25). Sim- 

344 ulated capture histories included 10 capture occasions, and data were generated by 

345 simulating true capture histories for individuals sequentially until 200 observed cap- 

346 ture histories were produced (each true history contributing either 0, 1, or 2 histories 
34? to the observed data). Survival and capture probabilities for each occasion were gen- 

348 erated from independent logit-normal distributions with means logit(.8) and logit(.5), 

349 and standard deviation .30. Recruitment rates were generated from independent log- 

350 normal distributions with mean log(.25) and standard deviation .30. The median 

351 number of true histories simulated before 200 observed histories were obtained was 

352 164 (min=150,max=180), the median number of unique individuals observed at least 

353 one time was 138 (min=127,max=148), and the median number of captures was 2 
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354 (min=l,max=10). 

355 Table [2] presents statistics comparing the mean-squared error (MSE) of the pos- 

356 terior means and the mean width and estimated coverage probability of the 95% 

357 credible intervals obtained from the alternative models. The MSE of the two-sided 

358 model and the combined-inference were similar for all parameters and smaller than 

359 those of the one-sided model by between 10% and 25%. Credible intervals for both 

360 the one-sided and two-sided model achieved the nominal coverage rate for all param- 

361 eters, but the credible intervals for the one-sided model were wider by approximately 

362 1 0%. In comparison, the credible intervals from the combined inference were narrower 

363 than those of the two-sided model by 20% or more but failed to achieve the nominal 

364 coverage rate. 

365 In the second scenario, we simulated data under the same model except that 

366 the capture probabilities were generated from a logit-normal distribution with mean 
36? logit(.5). The result is that capture probabilities were lowered and the number of 

368 observations per individual decreased. From the 100 data sets generated, the median 

369 number of true histories simulated before 200 observed histories were obtained was 196 

370 (min=170,max=235), the median number of unique individuals observed at least one 

371 time was 131 (min=119,max=143), and the median number of captures per observed 

372 individual was 1 (min=l,max=9). In this scenario, the MSEs of the posterior mean 

373 survival probabilities was similar between the two-sided model and the combined 

374 inference, but the two-sided model produced higher MSEs for both the recruitment 

375 and population growth rates. Credible intervals from the two-sided model were again 

376 narrower than those of the one-sided model by approximately 10%, but coverage 

377 probabilities from both models were above the nominal rate. Finally, although the 

378 credible intervals from the combined inference were narrower than those of the two- 

379 sided model by approximately 20%, the coverage of these intervals was also increased 

17 



380 and was near or above the nominal value. This result was due to the combination 

381 of non-informative priors and the relatively small amount of information in the data 

382 when capture probabilities are low; we discuss this point in Section [6] 

383 In the third scenario, we simulated data from the same model used in the first 

384 scenario except that both marks were seen with probability one each time an indi- 

385 vidual was captured (p^ = 1). This represents the extreme situation in which there is 

386 complete dependence between the two marks and no uncertainty in the true capture 
38? histories. Modeling the data from the left-side photographs, right-side photographs, 

388 or from both-sides using our model produces exactly the same results, and so we only 

389 compare the model of the left-side photographs with combined inference. The me- 

390 dian number of histories simulated in the 100 data sets before 200 observed histories 

391 were obtained was 215 (min=204,max=227) and the median number of captures per 

392 observed individual was 2 (min=l,max=10). In this scenario, the point estimates 

393 produced by the two models were almost exactly equal and the MSE of the two mod- 

394 els was indistinguishable. However, there were clear differences in the performance of 

395 the 95% credible intervals. While the intervals produced by combined-inference were, 

396 on average, 30% narrower the coverage of these intervals was well below the nominal 

397 value. 

398 [Table 2 about here.] 

399 5 Results 

400 The data provided in the ECOCEAN whale shark library contained a total of 96 

401 observed encounter histories for the 2008 study period. Of these, 27 histories (28%) 

402 were constructed from left-side photographs alone, 24 (25%) were constructed from 

403 right-side photographs alone, and 45 (47%) contained at least one encounter with 
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404 photographs taken from both sides simultaneously. Along with the model presented 

405 in Section [3j we computed inferences for p, /, and (p from the three alternative models 

406 described in Section HI 

407 Table [3] provides posterior summary statistics for the LBJS parameters model 

408 obtained from the two-sided model. Inferences about all parameters are relatively 

409 imprecise because of the relatively small number of individuals captured and the low 

410 capture probabilities, but the posterior means follow the expected patterns. Point 

411 estimates for the survival probability (the probability that a whale shark remains 

412 NMP between occasions) are at or above .90 in the first 2 periods, below .70 in the 

413 last two periods, and about .80 in between. The posterior mean recruitment rate is 

414 very high in week 2, suggesting that most of the sharks entered during this period, 

415 and lower thereafter. This table also provides summary statistics for the population 

416 growth rate, Xk = 4>k + fk, k = 1, . . . , K — 1, computed as a derived parameter. 

417 Although the 95% credible intervals for cover 1.00 for all k, the point estimates 

418 are greater than 1.00 for the first two periods, close to 1.00 in the next three periods, 

419 and less than .75 in the last two periods. This suggests that the aggregation of whale 

420 sharks grew during the first two periods, remained almost steady during the next 3 

421 periods, and declined during the last two periods. This supports the hypothesis that 

422 whale sharks aggregate at NMP to feed after the major coral spawn which occurred 



423 between April 9 and 12 in 2008 (Chalmers 2008, pg. 33). 



424 [Table 3 about here.] 

425 Table [4] provides posterior summary statistics for the conditional event probabili- 

426 ties. These results show that when sharks were encountered photographs were most 

427 often taken from both sides simultaneously (p% = .45(.36, .54)) and that the proba- 

428 bilities that an individual was photographed from either the left- or right-side only 
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429 were similar (p\ = .29(.20, .38) versus p\ = .21(.13, .29)). 

430 The posterior mean of n*, the number of unique sharks encountered during the 

431 2008 season, was 88 with 95% credible interval (82,93). The full posterior distribution 

432 of n* is shown in Figure [T] and compared with the prior distribution of n* generated by 

433 simulating data sets from the prior predictive distribution conditional on there being 

434 96 observed capture histories and at least 72 true histories (the minimum number 

435 given that 24 of 96 observed histories included right-side photographs alone). Whereas 

436 the prior distribution of n* is close to uniform the posterior distribution is strongly 

437 peaked and concentrates 95% of its mass between 82 and 93. We conclude that 

438 between 3 (3.1%) and 14 (14.6%) of the sharks encountered during the 2008 season 

439 were photographed from both the left- and right-sides on separate occasions without 

440 ever being matched. 

441 [Table 4 about here.] 

442 [Figure 1 about here.] 

443 Comparisons of the three chains starting from different initial values provided 

444 no evidence of convergence problems. Traceplots all indicated that the three chains 

445 converged within the burn-in period, GRB diagnostic values were all less than 1.02, 

446 and the estimated MCMC error was less than 2.6% of the posterior standard deviation 

447 for each parameter. Based on these results, we are confident that the chains converged 

448 and were run for sufficiently long to obtain reliable summary statistics. 

449 The plots in Figure [2] compare inferences for the survival, recruitment, and growth 

450 rates from the four alternative models. Posterior means from the four models are all 

451 very similar and the 95% credible intervals for all parameters overlap considerably. 

452 Comparison of the widths of the 95% credible intervals from the left- and right-side 
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453 data alone showed that the two-sided model provided improved inference in most, 

454 but not all, parameters. On average, the 95% credible intervals for the recruitment 

455 rates produced by the two-sided model were 93% and 69% as wide as those produced 

456 from the left- and right-side data alone. The 95% credible intervals for the survival 

457 probabilities produced by the two-sided model were 78% as wide as those from the 

458 right-side data, on average, but 103% as wide as those from the left-side data. This 

459 last result seems to be caused by issues with the upper bound on the survival prob- 

460 abilities as the 95% credible intervals for the logit transformed survival probabilities 

461 produced from the two-sided model were, on average, 90% and 89% as wide as those 

462 obtained from the left- and right-side data alone. Credible intervals produced via 

463 combined inference were on average 12% smaller than those obtained from the two- 

464 sided model; however, based on the results in the previous section, we believe that 

465 these intervals would not achieve the nominal coverage rate and do not reflect the 

466 variability of the parameters correctly. 

467 [Figure 2 about here.] 

468 6 Conclusion 

469 The simulation results presented in Section [4] illustrate the main advantages of our 

470 model over the previous approaches to analyzing mark-recapture data with multiple 

471 non-invasive marks In general, estimates from our model will be more precise than 

472 estimates based on only one mark in that posterior standard deviations will be smaller 

473 and credible intervals will be narrower but still achieve the nominal coverage rate. 

474 In contrast, the apparent gain from combining estimators computed separately for 

475 each mark under the assumption of independence is artificial and credible / confidence 

476 intervals computed by these methods will not achieve the nominal coverage rate - 
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477 particularly when the probability that both marks are seen simultaneously is high and 

478 the separate estimators are highly dependent. The differences between the methods 

479 are smaller when the capture probabilities are low (as in Simulation 2) but this is 

480 not surprising given that we have selected non-informative prior distributions. The 

481 posterior distribution is influenced more by the prior as the amount of information in 

482 the data decreases, and this would not have occurred if more informative priors had 

483 been selected. 

484 The disadvantage of the two-sided model is that the computations are more com- 

485 plex and more time consuming. The major challenge is that the number of possible 

486 true capture histories grows relative to the square of the number of observed capture 

487 histories. Possible solutions include deriving more efficient methods of computation 

488 or approximating the posterior distribution in some way. Further research is needed 

489 in this area before applying our model to large data sets. 

490 In the analysis of the whale shark data, we have focused on understanding the 

491 dynamics of the population (i.e., survival, recruitment, and population growths rates) 



492 but the method is easily adapted to estimate abundance. Following Link et al. (2010) 



493 one can include the null encounter history (vector of 0s) to the set of possible true 

494 histories and extend the vector of counts for the true histories to n* = (n*, . . . , n\, rig) 

495 where n* represents the number of individuals in the population that were never 



496 captured. The value n* = l T n* would then denote the total population size. Note 

497 that the prior on n* must also be modified by increasing M to allow n* to be larger 

498 than the number of observed histories. 

499 Our model can also be extended to combine data from any number of marks 

500 by expanding the set of constraints, and we expect that including more marks will 

501 strengthen differences between the alternative methods seen in the simulation study. 

502 The ratio between the median width of the credible intervals from the one-sided model 
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503 and the combined inference is close to 1.4 for all parameters in all three simulation 

504 scenarios - very near the a/2 reduction in standard deviation that occurs when two 

505 identically distributed and perfectly correlated random variables are averaged. Gen- 

506 erally, posterior standard deviations computed under the assumption of independence 
50? should decrease relative to the square root of the number of marks used, but the gain 

508 is artificial and the coverage of credible/confidence intervals will decrease further be- 

509 low the nominal value as more marks are used. In comparison, our model will account 

510 for the dependence between marks and will provide more precise inference without 

511 sacrificing coverage. 

512 One unusual aspect of the whale shark study is that individuals may be encoun- 

513 tered multiple times during each capture period. This would not occur if the capture 

514 periods were truly instantaneous, and in such cases both marks could only be observed 

515 in the same period if they were seen simultaneously. Event B could never occur. Our 

516 model could be fit to such data simply by specifying, a priori, that P(p4 = 0) = 1. 
si? and other researchers have begun to develop similar methods specific to such data 
sis (Conn and McClintock, pers. comm.). 
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Figure 1: Comparison of the prior and posterior distribution of n*. The prior dis- 
tribution of n* conditional on there being 170 observed capture histories is shown 
by the histogram with white bars. The posterior distribution of n* is shown by the 
histogram with grey bars. 
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Figure 2: Comparison between the two-sided model and the three alternative models. 
The plots on the left-side of the figure compare the posterior means (points) and 95% 
credible intervals (vertical lines) of the survival probability (top), recruitment rate 
(middle), and population growth rate (bottom) obtained from the four models. The 
plots on the right side of the figure display the posterior standard deviations from the 
three alternative models relative to the posterior standard deviation from the two- 
sided model. Results from the two-sided model are represented by the circles, from 
the left-side photographs only by the upward pointing triangle, from the right-side 
photographs only by the downward pointing triangles, and from combined inference 
by the diamonds. 



28 



Table 1: Example of possible observed and true capture histories. Suppose that 

the data comprises the six observed histories given in the top of the table. The 

possible true histories that may have generated this data include these six plus the 

four additional histories in the bottom of the table. 

Occasion 





1 


2 


3 


4 


5 


6 


7 


8 
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L 











2 














L 











3 








R 
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Table 2: Performance of the estimates from the three simulation studies. Each column 
of the table presents the MSE of the posterior mean relative to the MSE of the 
posterior mean of the one-sided model, and the median width and estimated coverage 
probability of the 95% credible intervals for the survival probability (0), recruitment 
rate (/), and growth rate (A) for one of the three models - one-sided (OS), two-sided 
(TS), or combined-inference (CI). Descriptions of the models are provided in the text 
in Section HI 

Simulation 1 Simulation 2 Simulation 3 

OS TS CI OS TS CI OS CI 



MSE 


1.00 


.89 


.87 


1.00 


.81 


.81 


1.00 


1.00 


Width 


.23 


.20 


.16 


.30 


.28 


.21 


.17 


.12 


Cover 


.97 


.96 


.90 


.98 


.99 


.93 


.95 


.84 


MSE 


1.00 


.88 


.81 


1.00 


.86 


.72 


1.00 


1.00 


Width 


.35 


.31 


.24 


.48 


.42 


.32 


.26 


.18 


Cover 


.97 


.95 


.90 


.99 


.98 


.96 


.95 


.84 


MSE 


1.00 


.88 


.82 


1.00 


.86 


.73 


1.00 


1.00 


Width 


.41 


.36 


.29 


.56 


.50 


.38 


.31 


.22 


Cover 


.98 


.97 


.95 


.99 


.99 


.97 


.97 


.87 
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Table 3: Posterior summary statistics for the demographic parameters <f>k, fk, and 
Pk obtained from the two-sided model. The columns of the table provide posterior 
means followed with equal-tailed 95% credible intervals. 

Occ (k) Survival ((pk) Recruitment (fk) Growth (A&) Capture (pk) 

1 0.90(0.67,1.00) 0.36(0.00,1.93) 1.26(0.76,2.83) 0.23(0.08,0.43) 

2 0.92(0.73,1.00) 2.40(0.08,6.41) 3.31(1.00,7.33) 0.19(0.05,0.33) 

3 0.82(0.54,1.00) 0.17(0.00,0.72) 0.99(0.64,1.56) 0.26(0.15,0.43) 

4 0.77(0.45,0.99) 0.09(0.00,0.36) 0.85(0.51,1.20) 0.22(0.13,0.34) 

5 0.82(0.49,1.00) 0.23(0.00,0.79) 1.05(0.63,1.65) 0.22(0.12,0.36) 

6 0.48(0.14,0.96) 0.06(0.00,0.29) 0.54(0.17,1.12) 0.25(0.14,0.42) 

7 0.66(0.16,0.99) 0.09(0.00,0.42) 0.75(0.20,1.28) 0.20(0.06,0.37) 

8 0.18(0.03,0.34) 
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Table 4: Posterior summary statistics for the conditional event probabilities. 

Event (j) Cond. Prob. (pj) 

1 0.29(0.20,0.38) 

2 0.21(0.13,0.29) 

3 0.45(0.36,0.54) 

4 0.06(0.01,0.13) 
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